e02ddf
e02ddf
© Numerical Algorithms Group, 2002.
Purpose
E02DDF Least-squares surface fit by bicubic splines with automatic knot
placement, scattered data
Synopsis
[lamda,mu,c,fp,rank,wrk,iwrk,nx,ny,ifail] = e02ddf(x,y,f,s,lamda,mu<,start,...
w,wrk,iwrk,nx,ny,lwrk,liwrk,ifail>)
Description
This routine determines a smooth bicubic spline approximation
s(x,y) to the set of data points (x ,y ,f ) with weights w , for
r r r r
r=1,2,...,m.
The approximation domain is considered to be the rectangle
[x ,x ]*[y ,y ], where x (y ) and x (y ) denote
min max min max min min max max
the lowest and highest data values of x (y).
The spline is given in the B-spline representation
n -4 n -4
x y
-- --
s(x,y)= > > c M (x)N (y), (1)
-- -- ij i j
i=1 j=1
where M (x) and N (y) denote normalised cubic B-splines, the
i j
former defined on the knots (lambda) to (lambda) and the
i i+4
latter on the knots (mu) to (mu) .
j j+4
The total numbers n and n of these knots and their values
x y
(lambda) ,...,(lambda) and (mu) ,...,(mu) are chosen
1 n 1 n
x y
automatically by the routine. The knots (lambda) ,...,
5
(lambda) and (mu) ,..., (mu) are the interior knots; they
n -4 5 n -4
x y
divide the approximation domain [x ,x ]*[y ,y ] into (
min max min max
n -7)*(n -7) subpanels [(lambda) ,(lambda) ]*[(mu) ,(mu) ],
x y i i+1 j j+1
for i=4,5,...,n -4; j=4,5,...,n -4. Then, much as in the curve
x y
case (see E02BEF), the coefficients c are determined as the
ij
solution of the following constrained minimization problem:
minimize
(eta), (2)
subject to the constraint
m
-- 2
(theta)= > (epsilon) <=S (3)
-- r
r=1
where: (eta) is a measure of the (lack of) smoothness of s(x,y).
Its value depends on the discontinuity jumps in
s(x,y) across the boundaries of the subpanels. It is
zero only when there are no discontinuities and is
positive otherwise, increasing with the size of the
jumps.
(epsilon) denotes the weighted residual w (f -s(x ,y )),
r r r r r
and S is a non-negative number to be specified by the user.
By means of the parameter S, 'the smoothing factor', the user
will then control the balance between smoothness and closeness of
fit, as measured by the sum of squares of residuals in (3). If S
is too large, the spline will be too smooth and signal will be
lost (underfit); if S is too small, the spline will pick up too
much noise (overfit). In the extreme cases the method would
return an interpolating spline ((theta)=0) if S were set to zero,
and returns the least-squares bicubic polynomial ((eta)=0) if S
is set very large. Experimenting with S-values between these two
extremes should result in a good compromise. Note however, that
this routine, unlike E02BEF and E02DCF, does not allow S to be
set exactly to zero: to compute an interpolant to scattered data,
E01SAF or E01SEF should be used.
Values of the computed spline can subsequently be computed by
calling E02DEF or E02DFF.
Parameters
e02ddf
Required Input Arguments:
x (:) real
y (:) real
f (:) real
s real
lamda (:) real
mu (:) real
Optional Input Arguments: <Default>
start (1) string 'c'
w (:) real ones(m,1)
wrk (lwrk) real if lower(start)=='c';wrk=...
... zeros(lwrk,1);end
iwrk (liwrk) integer if lower(start)=='c';iwrk=...
... zeros(liwrk,1);end
nx integer if lower(start)=='c';nx=0;end
ny integer if lower(start)=='c';ny=0;end
lwrk integer e02ddf18(length(x),...
... length(lamda),length(mu))
liwrk integer e02ddf20(length(x),...
... length(lamda),length(mu))
ifail integer -1
Output Arguments:
lamda (:) real
mu (:) real
c (:) real
fp real
rank integer
wrk (lwrk) real
iwrk (liwrk) integer
nx integer
ny integer
ifail integer